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We study the elasticity of fibrous materials composed of generalized stiff polymers. It is shown 
that in contrast to cellular foam-like structures affinc strain fields are generically unstable. Instead, 
a subtle interplay between the architecture of the network and the elastic properties of its building 
blocks leads to intriguing mechanical properties with intermediate asymptotic scaling regimes. We 
present exhaustive numerical studies based on a finite element method complemented by scaling 
arguments. 
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Cellular and fibrous materials are ubiquitous in na- 
ture and in many areas of technology. Examples range 
from solid or liquid foams over wood and bone to the 
protein fiber network of cells 0, S 01 • O n a mesoscopic 
level these materials are comprised of struts and mem- 
branes with anisotropic elastic properties. The systems 
differ widely in architecture. One finds patterns which 
are as regular as a honeycomb, as sophisticated as the 
particular design of a dragonfly's wing or simply ran- 
dom Q . The manifold combinations of architecture and 
elastic properties of the building blocks allows for a rich 
spectrum of macroscopic elastic responses. For cellular 
structures (see the foam- like structure in Fig. ^,) macro- 
scopic elasticity can already be understood by consider- 
ing the response of a single cell 0, ■ In these systems 
local stresses acting on an individual cell are the same as 
those applied on the macroscopic scale. In other words, 
the local deformation S of a cell with linear dimension l s 
follows the macroscopic strain 7 in an affine way such 
that it scales as 5 oc 7^ . This affinity of the deformations 
seems to be a reasonable approximation even for systems 
with some randomness in the size and the type of their 
unit cells [a, [3 • Since in affine models there can be no 
cooperativity between the elastic responses of individual 
cells, the effect of the assembled structure can simply be 
predicted by counting the numbers of cells. Fibrous net- 
works, on the other hand, are dramatically different al- 
ready in their morphology as can be inferred from Fig. ^ 
The presence of fibers introduces the additional meso- 
scopic scale of the fiber length I and generates a broad 
distribution of pore sizes that, in contrast to foams, has a 
non- vanishing weight even for the smallest cells 8] . This 
difference in architecture crucially affects the mechanical 
properties. Recently a non-affine regime has been iden- 
tified and characterized 0, ^| in two-dimensional 
networks of classical beams ("Mikado model") commonly 
used to model the mechanical properties of paper sheets 
[I2I Il3lll3. Il5| . The non- affinity of the deformation field 
necessarily implies that in these networks cooperativity 





FIG. 1: Illustration of the different architecture of (a) cellular 
and (b) fibrous materials. The foam in (a) is constructed by 
a Voronoi tessellation from the centers of the fibers in (b). 



effects play an important role. 

In this Letter we will address the question of how 
such non-trivial network geometry generates correlation- 
induced linear elastic behavior while cellular foam-like 
structures can generally be described by mean-field like 
approaches. By systematically tuning the force response 
properties of the individual elements we will be able to 
identify a new length-scale below which correlations drive 
the system away from the state of affine deformations. 
We will, moreover, describe the mechanism that gener- 
ates this length and calculate the resulting power law 
behavior of the elastic modulus by a scaling argument. 

The model is defined as follows. N anisotropic elas- 
tic elements, geometrically represented by straight lines 
of length I, are placed on a plane of area A — L 2 such 
that both position and orientation of the elements are 
uniformly random distributed. The length of the seg- 
ments, i.e. the distance l s between any two neighboring 
intersections, follows an exponential distribution || 



P(ls) = (Is) 
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with a mean value that is given in terms of the density 
p = Nl/A as (Z s ) = 7r/2p. At any intersection a perma- 
nent crosslink with zero extensibility is generated. This 
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constrains the relative translational motion of the two fil- 
aments, while leaving the rotational degrees of freedom 
independent (free hinges). Next to the architecture of the 
network we have to specify the linear elastic properties 
of the constituents. Previous studies 0, |n| have con- 
sidered classical beams of radius r and bending stiffness 
k. Loaded along their axis ("stretching") such slender 
rods have a rather high stiffness k s = 4n/l s r 2 , while they 
are much softer with respect to transverse deformations 
k±{l s ) — 3/t/Zf ("bending"). Here, we consider elastic 
elements where in addition to the mechanical stiffness of 
classical beams a more general stretching coefficient 
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is introduced. This may result from thermal fluctua- 
tions of the filament immersed in a heat bath of solvent 
molecules. The prefactor is chosen such that feii for a = 1 
reduces to the longitudinal entropic elasticity of a stiff 
polymer described by the wormlike chain model grafted 
at one end [l^. In this case the material length l p is 
called the persistence length of the polymer and quanti- 
fies the ratio of bending to thermal energy l p = n/ksT. 
The phenomenological exponent a allows us to extend 
our discussion to the broad class of systems for which 
fen is a monomial (with units energy per area) involving 
one additional material length l p . Having two longitudi- 
nal deformation modes the effective stretching stiffness is 
equivalent to a serial connection of the "springs" k s and 
ku. Setting l p — or, we can write k s cx fc|i(a = —2). The 
constant c is a material property of the specific polymer 
and has been chosen as c = 1.5 • 10 4 , which roughly cor- 
responds to the biopolymer F-actin. The precise value, 
however, is irrelevant what regards the thermal response 
fcii and only specifies the location of the cross-over to 
k s . The description of a thermally fluctuating network 
in terms of force constants ku , k s and k± is in the spirit 
of a Born-Oppenheimer approximation that neglects the 
fluctuations of the "slow variables", the cross-link posi- 
tions, while assuming the "fast" polymer degrees of free- 
dom to be equilibrated. By minimization of the internal 
energy with respect to the slow parameters, the modulus 
for a given macroscopic deformation can be calculated. 
In practice, we apply a shear strain of 7 = 0.01 and thus 
determine the shear modulus G using periodic boundary 
conditions on all four sides of the simulation box. The 
minimization procedure is performed with the commer- 
cially available finite element solver MS CM ARC. 

To set up a reference frame, it is instructive to discuss a 
mean-field approach which assumes that correlations be- 
tween neighboring segments are absent. As outlined in 
the introduction this corresponds to a foam-like behavior, 
where the response is fully described by the properties of 
an average segment of length (l s ) cx p^ 1 . Marking the 
force constants of this segment by an overbar, we can ex- 
press them in the form (neglecting numerical prefactors) 




FIG. 2: Shear modulus Gl 3 / K as function of persistence length 
lp/l for various densities pi. The second branch in the upper 
right corner (pi > 120) is obtained by suppressing the me- 
chanical response ("fc s — > 00"). The dashed line indicates the 
three regimes as obtained by Eq. 



k± ~ up 3 , fen ~ k±(pl p ) a and k s ~ k±(pr)~ 2 , respec- 
tively. The corresponding deformation modes will act 
as springs connected in series .5] such that the modulus 
takes the form 
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The network will thus show a crossover from thermal 
stretching to bending at l v ~ (l s ) and to mechanical 
stretching at r w (l s ). This behavior is indicated by 
the dashed line in Fig. [21 where it can be compared with 
the actual results of our numerical analysis. The nor- 
malized shear modulus Gl 3 /k is shown as a function of 
dimensionless persistence length l p /l for a set of dimen- 
sionless densities pi. At large l p /l (in the right part of 
the plot) we recover purely mechanical behavior charac- 
terized by G oc fc.s consistent with the mean-field picture 
of Eq. (|3J) [jJJ llll ll2| . Our main interest, however, lies 
in the regime of small l p /l (in the left part of the plot), 
where the persistence length is small enough for thermal 
fluctuations to become relevant. To analyze the modulus 
in the thermal regime (where k s is irrelevant) it will be 
helpful to use dimensional analysis and write the mod- 
ulus in terms of the two remaining response coefficients 
k± and feii of an average segment 



G(k, I, l pi p) = k±g(pl, k\\/k±) . 



(4) 



The first argument of the scaling function g, the density 
x = pi, is of geometrical origin and counts the num- 
ber of crosslinks per filament. The second argument, 
y = k\\/k± ~ pl p , relates to the energy balance between 
stretching and bending of an average segment and marks 
a crossover at j/ » 1 or l p « (l s ). From Fig. and the 
inset of Fig.|3one infers that for low densities g — yf(x), 
implying for the modulus G = k\\f(pl). This linear de- 
pendence on the "pre-averaged" stretching compliance 
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FIG. 3: Scaling function g as a function of pl p for various 
x — pi > 60; Inset: The same for x < 60. 



hints at a foam-like stretching dominated regime |l7j | 
where correlations are absent. As one can also infer from 
these figures the domain of validity of this linear regime 
is extremely narrow and confined to low densities x < 20 
and persistence lengths y <C 1 . For medium and high den- 
sities Fig. |21 shows two non-trivial scaling regimes where 
g(x 1, y) = y z becomes independent of x (and there- 
fore of the filament length I) and exhibits power law be- 
havior with exponents z = 0.46 and z = 0.9 for small and 
large values of y, respectively. In both cases, the modu- 
lus can be written in the form of a generalized geometric 
average 



G oc k]~ z k 



(5) 



which has to be contrasted with Eq. @, where bending 
and stretching modes are assumed to superimpose lin- 
early. Here, correlation effects between the segments in- 
duce the non-trivial form of the modulus and distinguish 
our fiber network from the ordinary foam-like behav- 
ior obtained by single segment considerations. Whereas 
foams may be considered as a limit where the filament 
length is zero (or rather identical to the cell size) the scal- 
ing limit of fiber networks corresponds to infinite fiber 
length. 

To understand the origin of the correlations one has to 
take into account the full distribution of segment lengths, 
Eq. QJ. This will have a pronounced effect on an affine 
deformation field <5 a g oc jl s as can be seen by consider- 
ing the axial force /|| along an affinely stretched segment 
of length l s , /|| = fc||6 a ff ~ nlp^/l^ ". In any but the 
purely mechanical situation, where a = — 2 (and thus 
f s ~ K'f/r 2 ), /|| strongly depends on the segment length. 
This implies that, in general, two neighboring segments 
on the same filament produce a net force at their com- 
mon node that has to be taken up by the crossing filament 
which then, preferentially, will start to bend. From the 
exponential distribution of segment lengths in Eq. one 



can easily show that the size of these residual forces Sf 
can be arbitrarily large. The corresponding probability 
distribution Q(Sf) shows polynomial (fat) tails 



Q(Sf)^5f-^ 3 P(l s ),6f 
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and has a diverging mean value since P(l s = 0) ^ 0. 
As a consequence there are always residual forces high 
enough to cause significant bending of the crossing fila- 
ment. Hence we conclude that an affine deformation field 
is unstable and that the system can easily lower its energy 
by redistributing the stresses to relieve shorter segments 
and remove the tails of the residual force distribution 

QW). 

This mechanism can be used to derive an expression 
for the modulus in the parameter region y -C 1, where 
the value of the exponent z — 0.46 indicates that bend- 
ing and stretching deformations contribute equally to the 
elastic energy. We assume that segments up to a critical 
length l c - to be determined self-consistently - will fully 
relax from their affine reference state to give all their en- 
ergy to the neighboring segment on the crossing filament. 
The energy of segments with l s > l c will then have two 
contributions. First, a stretching part from the imposed 
affine strain field (for simplicity, we will set a = 1 in what 
follows) 
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Second, a bending part 



w b (l s ) ~ k±S 2 s ~ /«7 
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that only arises if the segment under consideration is 
neighbor to an element on the crossing filament with 
l' s < l c (the prime refers to the neighboring small seg- 
ment). Adding both contributions and averaging over 
all segments / s > l c and l' s < l c we arrive at the ex- 
pression w ~ kj 2 p(pl p /x c + x c ), where x c '■= pl c <C 1 in 
the parameter range of interest. Minimizing with respect 
to x c gives the required expressions x™ ln ~ (Ipp) 1 ^ 2 and 
G ~ p 2 w min / r y 2 ~ np 7 / 2 ^ 2 corresponding to a value 
z = 1/2 for the exponent that compares well with the 
measured value z = 0.46. Repeating the calculation for 
general values of a gives z(a) = a/(l + a). We have 
verified this result by simulations with an accuracy of 
about ten per cent 181. The non-trivial behavior of G 
observed in Figs. El and can thus be explained by a 
length scale l c = (l s )(l p / {Is)) 1 / 2 below which the affinity 
of the deformation field breaks down. The mechanism is 
illustrated in Fig. where a histogram for the fraction 
of energy stored in segments of various lengths is shown. 
Increasing the persistence length, the short segments one 
after the other loose their energies in favor of additional 
excitations in longer segments. 
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FIG. 4: Fraction of energy stored in the various segment 
lengths; the curves correspond to different persistence lengths 
at a density of pi — 80, equivalent to (l s )/l ~ 2 • 10~ 2 . 



mixed regime we found a stretching dominated regime 
that is characterized by exponents similar to those ob- 
tained in the foam-like stretching regime. Microscopi- 
cally, however, the stress is distributed in a completely 
different way and is dominated by very few, very long 
fiber segments. As mentioned above, one particular ap- 
plication of our model is to the macroscopic linear re- 
sponse of stiff polymer networks. The results may there- 
fore be directly relevant for two-dimensional networks of 
the filamentous biopolymer F-actin, assembled on top 
of micro-fabricated pillars [llSj . In addition, it might 
shed new light on very recent rhe olog ical measurements 
on cross- linked actin networks |2(i l21| , which emphasize 
the single-polymer origin of the measured elastic moduli. 
Our simulations, on the contrary, highlight the poten- 
tially non-trivial effects of inter-polymer correlations on 
the macroscopic elasticity. 



When, eventually, l c w l p w (J s ) (y « 1) the affine 
strain field does not serve as a reference configuration 
any more, since it is strongly perturbed by a majority of 
segments with l s < l c . Moreover, the unloading of the 
smaller segments now produces significant stretching de- 
formations of their neighbors on the same filament such 
that the available energy for bending of the crossing fil- 
ament is reduced. At this stage, one enters the second 
intermediate asymptotic regime where, like in the affine 
regime at low densities, stretching modes dominate the 
modulus. As can be seen in Fig. 21 only the longest 
segments carry substantial amounts of energy such that 
the displacement field must be highly non-affine. For 
lp/l > 1.5 10~ 2 about 90% of the energy is stored in the 
longest 30% of the segments. In this parameter range, 
bending is on average the softer mode y = k\\/k± 3> 1 
and therefore contributes only very little to the total 
energy. Raising the density to still higher values, it is 
conceivable from our data that the exponent z — 0.9 ap- 
proaches z — 1, which would mean that the energy in the 
bending modes is completely negligible and G ~ fc|| as in 
the affine regime. A transition into a regime dominated 
by the low-energy bending modes would not be favor- 
able, however. As is known from the mechanical fiber 
model jlflj |. such a regime must not be characterized by 
the pre-averaged force constant k± = k±((l s )) but by an 
effective stiffness (kj_) ~ k/^ 3 with a new length scale 
£ = l(pl)~^l z and /i «s 6.7 that is highly dependent on 
fiber length 

In summary, we found that for a broad range of pa- 
rameters the macroscopic shear modulus is asymptoti- 
cally independent of the fiber length. Affine stretching 
is energetically unstable towards a redistribution of en- 
ergies in favor of longer segments. This gives rise to a 
correlation-induced elasticity that cannot be explained 
within a "single cell" picture. In addition to a non-trivial 
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